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The spatial structure of populations is a key element in the understanding of the large scale spreading 
of epidemics. Motivated by the recent empirical evidence on the heterogeneous properties of transportation 
and commuting patterns among urban areas, we present a thorough analysis of the behavior of infectious 
diseases in metapopulation models characterized by heterogeneous connectivity and mobility patterns. We 
derive the basic reaction-diffusion equation describing the metapopulation system at the mechanistic level 
and derive an early stage dynamics approximation for the subpopulation invasion dynamics. The analytical 
description uses degree block variables that allows us to take into account arbitrary degree distribution of 
the metapopulation network. We show that along with the usual single population epidemic threshold the 
metapopulation network exhibits a global threshold for the subpopulation invasion. We find an explicit 
analytic expression for the invasion threshold that determines the minimum number of individuals traveling 
among subpopulations in order to have the infection of a macroscopic number of subpopulations. The 
invasion threshold is a function of factors such as the basic reproductive number, the infectious period 
and the mobility process and it is found to decrease for increasing network heterogeneity. We provide 
extensive mechanistic numerical Monte Carlo simulations that recover the analytical finding in a wide 
range of metapopulation network connectivity patterns. The results can be useful in the understanding of 
recent data driven computational approaches to disease spreading in large transportation networks and the 
effect of containment measures such as travel restrictions. 
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1. Introduction 

The metapopulation modeling approach is an essen- 
tial theoretical framework used in population ecology, ge- 
netics and adaptive evolution to describe population dy- 
namics whenever the spatial structure of populations is 
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known to play a key role in the system's evolution (Han 
ski Gilpinl [19971 [Hanski GaggiottH [20041 [Tilman 



Kareiva 1997 



Bascompte fc Sole[|1998 ). Metapopulation 



ley, 2007). The arrival of the infection in any subpopu- 
lation and its epidemic evolution are determined by the 
coupling generated by the mobility processes among sub- 
populations. The metapopulation dynamics of infectious 
diseases has generated a wealth of models and results 
considering both mechanistic approaches taking explic- 



models rely on the basic assumption that the system un- 
der study is characterized by a highly fragmented environ- 
ment in which the population is structured and localized 
in relatively isolated discrete patches or subpopulations 
connected by some degree of migration. Classic metapop- 
ulation dynamics focuses on the processes of local ex- 
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tive coupling approaches where the diffusion process is 
expressed as a force of infection coupling different sub- 



tinction, recolonization and regional persistence (Levins 



1969 1970), as the outcome of the interplay between mi- 
gration processes among unstable local populations and 
population dynamics (e.g birth and death rates, competi- 
tion and predations). This paradigm is extremely useful 
also in the case of infectious diseases, and can be ap- 
plied to understand the epidemic dynamics of spatially 
structured populations with well defined social units (e.g. 
families, villages, city locations, towns, cities, regions) 



populations ( Bolker fc Grenfell 


1995 


Lloyd fc Mayl 


19961 


Earn et aL\ 19981 Rohani et al. 


1999 


KeelingI [20001 [Park[ 


et al\ 2002 Vazquezl 2007|. 


Recently, the metapopu- 



lation approach is being revamped in computational ap- 
proaches for the large scale forecast of infectious disease 
spreading (Grai s a/.] |2004 ; "Hufna gel efal\ [2004[ [Col- 
izza et al., 2006a]' Cooper et a/. , 2006^|Colizza et 



Hollingsworth et al.^ 2006 Riley 2007) 



2007^ 



connected through individual mobility (Hethcote 



May fc Anderson] [19791 [Anderson fc May[ |1984| [May fc l 



1978 



Ander son', '1984'; 'Bolke r fc Grenfell[ [19931 11995| [K eeling ' 
Rohani, 2002, Lloyd fclTay] [1996] [Grenfeh fc Harwood, 



Metapopulation epidemic models, especially at the 
mechanistic level, are based on the spatial structure of the 
environment, and the detailed knowledge of transporta- 
tion infrastructures and movement patterns. The increas- 
ing computational power and informatics advances are 
beginning to lift the constraints limiting the collection of 
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large spatiotemporal data on human behavior and demog- 
raphy, finahy ahowing for the formulation of realistic data 
driven models. On the other hand, networks which trace 
the activities and interactions of individuals, social pat- 
terns, transportation fluxes, and population movements 
on a local and global scale ( Liljeros et al] 
berger et al] [20041 [Barrat et aL\ |2QQ4| |Guimera et aL\ 



2005 



Chowell et al 2003 ) have been analyzed and found 
to exhibit complex features encoded in large scale het- 
erogeneity, self-organization and other properties typical 



of complex systems ( 


Albert & Barabasi, 2002; Dorogovt- 


sev & Mendesl 2003 


Newmanl |2003t Pastor-Satorras & 


Vespignani| 2004). In particular, a wide range of societal 



and technological networks exhibits very heterogeneous 
topologies. The airport network among cities (Barrat 



et a/., 2004 Guimera et al. 2005), the commuting pat- 



terns in inter and intra- urban areas ( Chowell et al. 2003 



Barrett et al] |2000l |De Montis et al.\ |2007D, and sev- 



eral info-structures rtPastor-Satorras & Vespignani 2004) 



are indeed characterized by networks whose nodes, rep- 
resenting the elements of the system, have a wildly vary- 
ing degree, i.e. the number of connections to other ele- 
ments. These topological fluctuations are mathematically 
encoded in a heavy-tailed degree distribution de- 
fined as the probability that any given node has degree 
/c, and have been found to have a large impact on epi- 



demic phenomena on complex contact patterns (Ander- 
son &: May[ [19921 [Pastor-Satorras &: Vespignani] |2001a|6| 



et 



Moreno et 'al] [2002] [Lloyd Ma^ [2501] [Barthelemyl 



2005[ ). 

Motivated by the above findings we provide here the 
analysis of the behavior of epidemic models in metapopu- 
lation systems with heterogeneous connectivity patterns. 
In order to have a mechanistic description of the system, 
we derive the deterministic reaction-diffusion equations 
describing the evolution of the epidemic in the metapopu- 
lation systems. The heterogeneity of the network is taken 
explicitly into account by introducing degree block vari- 
ables that provide results expressed as functions of the 
moments of the degree distribution of the substrate net- 
works. In order to account for the discreteness of the 
system and microscopic fluctuations in the diffusion pro- 
cesses we derive also coarse grained equations for the inva- 
sion dynamics at the subpopulation level. The system is 
characterized by the standard (i.e. single population) epi- 
demic threshold and by a global invasion threshold pro- 
viding the condition for the infection of a macroscopic 
number of subpopulations. The first threshold defines the 
usual reproductive number Rq > 1 that is just a function 
of the disease parameters while the second threshold de- 
fines a subpopulations reproductive number > 1 that 
depends also on the diffusion rate of individuals among 



subpopulations i\Bal\\ [1997[ [Cross et al] [2005[ [2007] ) . We 
find an explicit analytic expression in the limit Rq ^ 1 
for the invasion threshold that is found to depend also on 
the network heterogeneity. The larger is the network het- 
erogeneity and the smaller is the diffusion rate that guar- 
antees the invasion of a finite fraction of subpopulations. 



This result provides a framework for the understanding of 
the evidence collected on the interplay between travel and 



global spread of infectious diseases (Viboud et al. 2006) 
and the poor effectiveness of travel restrictions in the con- 



tainment of epidemics ( 


Cooper et 


al 


. 2006 


Hollingsworth 


et al. , 2006 


Colizza et al. , 2007a 


). Finally, the analytic 



results are confirmed by mechanistic Monte Carlo simu- 
lations for the infection dynamics in the metapopulation 
system, in which each single individual is tracked in time 
to account for the discreteness of the processes involved. 
Heterogeneous connectivity patterns among subpopula- 
tions are modeled and different values of the parameters 
involved are considered to validate the theoretical results. 

The paper is organized as follows. Section 2 introduces 
the metapopulation epidemic model on a heterogeneous 
network of connections among subpopulations. Two dif- 
ferent kinds of mobility processes are introduced in Sec- 
tion 3 to analyze the stationary diffusion properties of 
the system. Section 4 incorporates the mobility processes 
analyzed in the previous section into a metapopulation 
epidemic model. Stochastic effects and discrete descrip- 
tion of the processes are considered with a tree-like ap- 
proximation for the analysis of the invasion dynamics at 
the level of the subpopulations. The effect of diffusion 
properties on the invasion dynamics are analyzed and re- 
lated to the existence of an invasion epidemic threshold 
for the metapopulation system. In Section 5 the behav- 
ior of the system above the invasion threshold is studied 
by mechanistic reaction-diffusion equations using a deter- 
ministic degree block variables representation. Finally, in 
Section 6 we report extensive mechanistic Monte Carlo 
simulations which confirm the analytical findings of the 
previous sections. 



2. Metapopulation mechanistic model as a 
microscopic reaction-diffusion process 

Metapopulation models describe spatially structured 
interacting subpopulations, such as city locations, urban 



areas, or defined geographical regions (Hanski & Gag- 



giotti, 2004 Grenfeh & Harwood, 1997). Individuals 



within each subpopulation are divided into classes denot- 



ing their state with respect to the modeled disease (An 



derson & May 1992) — such as infected, susceptible, im- 
mune, etc. — and the compartment dynamics accounts for 
the possibility that individuals in the same location may 
get into contact and change their state according to the 
infection dynamics. The interaction among subpopula- 
tions is the result of the movement of individuals from one 
subpopulation to the other. It is clear that the key issue 
in such a modeling approach is how accurately we can 
describe the commuting patterns or traveling of people. 
In many instances even complicate mechanistic patterns 
can be accounted for by effective couplings expressed as a 
force of infection generated by the infectious individuals 
in subpopulation j on the individuals in subpopulation 
i ([Bolker & Grenfell| [19951 [Lloyd & May| [19961 [Earn et 



alj [19981 [Rohani et aL\ [19991 [Keelmgl [2QQQ1 [Fark aL\ 
2002). More realistic descriptions are provided by ex- 



plicit mechanistic approaches which include the detailed 
rate of traveling/commuting obtained from data or from 
empirical fit to gravity law models (for a recent reference, 
Viboud et al ( 2006[ )), accompanied by the associated 



see 



mixing subpopulations Nij denoting the number of indi- 
viduals of the subpopulation i present in the subpopu- 



1995). 



lation j (Keeling & Rohani 2002 Sattenspiel & Dietz 



A simplified mechanistic approach uses a markovian 
assumption in which at each time step the movement of 
individuals is given according to a matrix dij that ex- 
presses the probability that an individual in the subpop- 
ulation i is traveling to the subpopulation j. The marko- 
vian character is in the fact that we do not label individ- 
uals according to their original subpopulation (e.g. home 
in a commuting pattern framework) and at each time step 
the same traveling probability applies to all individuals in 
the subpopulation without having memory of their origin. 
This approach is extensively used for very large popula- 
tions in the case the traffic Wij among subpopulations 
is known by stating that dij ~ Wij/Nj. Several model- 
ing approaches to the large scale spreading of infectious 
disease (Baroyan et a/.| |1969| [Rvachev & Longini, 1985 
Longinil |1988| [Fla hault & Valler on| [19911 |Gr ais eFot 



2om 



20041 [Hufnagel et a/.H2004, [Colizza et a/.||2006a|6 



2007a) use this mobility process based on transportation 
networks for which it is now possible to obtain detailed 
data. 

In their simplest formulation markovian mechanistic 
mobility processes are equivalent to the classic reaction 
diffusion processes used in many physical, chemical and 



pen 



biological processes (Marro & Dickman, 1999 van Kam- 



1981 Murray 2005). The reaction-diffusion frame- 



work (Colizza et a/., 20076) considers that the occupa- 



tion numbers A^^ of each subpopulation can have any in- 
teger value, including Ni = 0, that is, void nodes with 
no individuals. The total population of the metapopu- 
lation system is N = A^^ and each individual diffuse 
along the edges with a diffusion coefficient d^j that de- 
pends on the node degree, subpopulation size and/or the 
mobility matrix. A sketch of the metapopulation model 
which shows the different scales of the system is shown 
in Fig. [T] The system is composed of a network sub- 
strate connecting subpopulations over which individuals 
diffuses. Each subpopulation is represented by a node i of 
the network. We consider that each node i is connected 
to other ki nodes according to its degree resulting in a 
network with degree distribution P{k) and distribution 
moments (/c") = k'^P{k). 

In the case of large metapopulation systems with a 
high level of heterogeneity the analytical description of 
the metapopulation model in terms of specific features of 
each single subpopulation is extremely complicate. In the 
following we propose an analytical framework that uses 



degree block variables to obtain the dynamical equations 
describing the system's behavior, relying on the empirical 
evidence pointing to a statistical equivalence of subpopu- 
lations having the same degree. 
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Figure 1: Schematic representation of a metapopulation 
model. The system is composed of a heterogeneous network of 
subpopulations or patches, connected by migration processes. 
Each patch contains a population of individuals who are char- 
acterized with respect to their stage of the disease (e.g. suscep- 
tible, infected, removed), and identified with a different color 
in the picture. Individuals can move from a subpopulation to 
another on the network of connections among subpopulations. 



2.1. METAPOPULATION NETWORKS WITH 
HETEROGENEOUS TOPOLOGY 

In the real world, the network specifying the cou- 
pling between different subpopulations is in many cases 
very heterogeneous. Examples can be drawn from several 
transportation infrastructures, commuting data and cen- 



sus information (Chowell et al\ 2003 Barrett et al\ 2000 



Barrat et a/.] |2004| [Guimera a/.|[2005|[De Montis al 
2007*). A particularly relevant one in the field of epidemic 
modeling is given by the airline transportation network. 
In this case the coupling is provided by the number of 
passengers traveling on a given route connecting two air- 
ports, thus yielding a transfer of individuals between the 



corresponding urban areas. For instance, [Barrat et al. 
(2004) reports a detailed study of the International Air 



Transport AssociatiorQ^^^^base which contains the com- 
plete list of world commercial airport pairs connected by 
direct fiights. Moreover, to each direct fiight connection 
between airports j and ^ is assigned a weight Wj^ which 
corresponds to the number of available seats or passengers 
on the given route. The obtained network displays high 
levels of heterogeneity both in the connectivity pattern 
and in the traffic capacities, as revealed by the broad dis- 
tributions of the number of connections of each airport, 
of the travel fiows between connected airports and of the 
traffic in terms of number of passengers handled by each 
airport ( Barrat et al. 2004 ) . These results have been con- 



firmed on the US subnetwork along different years and 
considering both market and segment traffic data^ and 
analogous results are recovered by analyzing commuting 
patterns data, intra-city traffic among locations, and sev- 
eral other data sets concerning the movements of people 



*IATA, International Air Transport Association, http: / / www.iata.org/| 
^BTS, Bureau of Transportation Statistics, http: / /www. bts.org/ 
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and goods dChowell et al] |2QQ3| [Barrett et fl/.[|2QQQ[ [Bar- 
rat et g/.l [2QQ4JjGuimera a/.[ |2QQ5| |De Montis et aL \ 



2007). In Fig. ^ we report the degree and weight prob- 
abihty distributions in some examples of these networks. 
In many cases we find heavy-tailed distributions varying 
over several orders of magnitude. For instance the airline 
traffic among different urban areas in the world shows a 
probability distribution P{w) - where w is the traffic on 
a single connection - varying over six orders of magnitude 
(see Fig. [2I 
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Figure 2: Degree (left column) and weight (right column) 
probability distributions for three different datasets of the 
movement of people and goods: (A)-(B) world-wide airport 
network, where the weight represents the number of seats on 
the flights connecting two different airports (source: lATA); 
(C)-(D) US Air Domestic Market, where the weight represents 
the number of passengers flying on a given origin-destination 
itinerary (source: BTS); (E)-(F) US Air Domestic Market for 
freight transportation, where the weight represents the amount 
of freight (expressed in lb) transported from an origin airport 
to its flnal destination airport (source: BTS). All datasets 
show heavy-tailed distributions both in the number of connec- 
tions and in the amount of people/good transported. 

In addition, it is also possible to find some general 
statistical laws relating the traffic and the degree of each 
node in the network. In general the behavior of the aver- 
age weight along the connection between two subpopula- 
tions with degree k and k' is a function of their degree 



{wkk') = wo{kk'f 



(1) 



where wq and the exponent depend on the speciffc sys- 
tem (e.g. ^ :^ 0.5 in the world-wide air transportation 



network (Barrat et al. 2004)). A related quantity is the 



total average traffic of the subpopulations with degree 
k that behaves as 



(2) 



Here the proportionality constant A and the exponent 
1 + ^ are deffned by the sum rule T/^ = WQ{kk'Y that 



must be satisffed on average. This last relation gauges 
the proportionality constant yielding A — {k^^^)wQ/ {k) . 
It is important to stress that the above relations deffnes a 
statistical equivalence of the subpopulations of degree k. 
In the following we want to address the questions of how 
the large scale complex features (scale invariance, extreme 
heterogeneity, unbounded ffuctuations) of interaction and 
communication networks affect the behavior of metapop- 
ulation models by deffning a mechanistic approach based 
on the master equations describing the disease dynamics 
as a microscopic reaction-diffusion process. The equiva- 
lence assumption that will be used throghout the rest of 
this paper is crucial in order to carry out the analytical 
treatment of the model. While this assumption is indeed 
recovered in several data sets, exceptions and ffuctuations 
have been noted that would require more complicate cal- 
culation schemes. 



3. Mobility processes and diffusion properties in 
heterogeneous networks 

In order to tackle the description of metapopulation 
models at the mechanistic level, let us ffrst analyze the 
simple diffusion process of a global population of N indi- 
viduals who diffuse in a network made of V nodes, each 
representing a subpopulation. Each node i stores a num- 
ber A^^ of individuals as deffned in section 2. In order to 
take into account the topological ffuctuations of the net- 
work we have to explicitly consider the presence of nodes 
with a widely ffuctuating degree k. A more convenient 
representation of the system is therefore provided by the 
quantities 

= ^ E ' 



where Vk is the number of nodes with degree k and the 
sums run over all nodes i having degree ki equal to k. The 
variable is therefore representing the average number 
of individuals in subpopulations within the degree block 
k. This representation is assuming that all subpopula- 
tions of the same degree are statistically equivalent. 

Let us assume a general framework in which the in- 
dividuals move from a subpopulation with degree k to 
another with degree k' with a diffusion rate dkk' that de- 
pends on the degrees of the origin and destination subpop- 
ulations. The probability of leaving a subpopulation with 
degree k is then given by pk = XI /c' ^(^1^)^^^'' where 
P{k'\k) is the conditional probability that any given edge 
departing from a node of degree k is pointing to a node 
of degree k' . In the following, we will ffrst write the 
equations for the dynamics of the individuals under this 
generic type of diffusion, and then address two speciffc 
diffusion rates to ffnd the stationary solutions. 

The dynamics of individuals is simply represented by 
a mean-ffeld dynamical equation expressing the variation 
in time of the subpopulations (t) in each degree block. 
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This can be easily written as: 

dtNkit) = -pkNk(t) ^kY,nk'\k)dk'kNk'(t). (4) 

k' 

The first rhs term of the equation just considers that only 
a fraction of particles pk moves out of the node. The sec- 
ond term instead accounts for the particles diffusing from 
the neighbors into the node of degree k. This term is 
proportional to the number of links k times the average 
number of particles coming from each neighbors. This is 
equal to average over all possible degrees k' the fraction 
of particles moving on that edge dk'kNk'{t) according to 
the conditional probability P{k'\k). In the following we 
will consider the case of uncorrelated networks in which 
the conditional probability does not depend on the orig- 



inating node, i.e. P{k'\k) = k'P{k')/{k) ( 


Dorogovtsev & 


Mendes 2003 Pastor-Satorras a/.l|2001 


) . This relation 



simply states that any edge has a probability to point to 
a node with degree k' that is proportional to the degree of 
the node. By using this form for P{k'\k)^ the dynamical 
rate equation Q for the subpopulation densities reads as 



dtNk(t) = -pkNk(t) ■ 



{k) 



Y,k'P{k')dk'kNk'{t). (5) 



In the following subsections we solve the previous set of 
equations for different diffusion processes that consider 
diffusion rates depending on the traffic of each node or on 
the population size of each subpopulation. 



3.1. TRAFFIC DEPENDENT MOBILITY RATES 

Here we assume that the probability an individual 
leaves a given subpopulation is independent of its degree 
k^ Pk = p \/k. If we also assume homogeneous diffusion 
along any given connection, individuals have the same 
probability to move along anyone of the links departing 
from the node at which they are located. In this case the 
diffusion rate along any given link of a node with degree 
k will be simply equal to 



dkk' = p/k. 



(6) 



This is obviously not the case in a wide range of real 
systems where the extreme heterogeneity of traffic is well 
documented (see subsection 2.1). A more realistic process 
therefore considers the movement of individuals to be pro- 
portional to the traffic intensity along a given edge. This 
is simply obtained by defining a heterogeneous diffusion 
probability for any given individual to go from a subpop- 
ulation of degree /c to a subpopulation of degree k' as 



dkk' P 



woikk'Y 



(7) 



This relation states that the diffusion rate p is still con- 
stant in each subpopulation but the individuals move on 
each connection in a proportion dependent on the actual 



traffic on the connection. The denominator = 
provides the correct normalization in order to ensure that 
by summing over all k edges departing from the node the 
overall diffusion rate is p. 

By using the expression of eq. ^ for dkk' and im- 
posing pk = p, the dynamical rate equation Q for the 
subpopulation densities reads as 



dtNk{t) = -pNk{t)^pk^'^'^ 



A{k)^ 



P{k')Nu'(t). (8) 



The stationary solution dtNj.{t) = does not depend 
upon the diffusion rate p that just fixes the time scale 
at which the equilibrium is reached and has the solution 

Wo 



A{k) 



(9) 



where TV = P{k)Nk{t) represents the average sub- 
population size. The explicit form of the normalization 
constant A = {k^~^^)wo/ {k), finally provides the explicit 
stationary solution 

The above solution states that the population of each 
node scales with the node degree in the stationary limit. 
The above behavior is simply the effect of the diffusion 
process that brings a large number of individuals in well 
connected, high traffic nodes, thus showing the impact of 
network's topological (i.e. dependence on k) and traffic 
(i.e. dependence on 6) fluctuations on the individuals den- 
sity behavior. When ^ = we recover the homogeneous 
diffusion case in which d^k' = dk = p/k^ obtaining 

= A«, (II) 

In this case the subpopulation density is just fixed from 
topological fluctuations and the exponent 6 clearly ap- 
pears as the parameter that takes into account the traffic 
fluctuations. It is worth remarking that in this frame- 
work, the subpopulation size as a function of the degree 
is constrained by the diffusion processes, a feature that 
has not to be expected in real systems where the popula- 
tion size of local patches can be considered as an indepen- 
dent variable. On the other hand, the degree dependence 
is close to those observed in real systems where in sev- 
eral cases it is possible to find a relation ~ k^ with 
0.5 < 



< 1.5 (Colizza et 



2006 a^6). 



We can thus 



consider the obtained stationary state as a first approxi- 
mation to the real case and use the exponent 6 to explore 
different levels of heterogeneity. 

3.2. POPULATION DEPENDENT MOBILITY RATES 

In a more general perspective, it is important to have 
the possibility of considering the population densities Nj. 
as independent variables. This is indeed the case of many 
metapopulation models in which the diffusion process rep- 
resents the travel of individuals between subpopulations. 
In this framework the number of people traveling from a 
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subpopulation to the other in a unitary time scale is a de- 
fined number Wij and the number of travehng individuals 
is independent from the population size A^^. This amounts 
to state that each individual in the subpopulation has a 
diffusion rate Wij/Ni where Wij is the total num- 
ber of people traveling out of city i in the unitary time 
scale. In other words, the diffusion rate of each individual 
is inversely proportional to the population size. In order 
to have a non-pathological stationary state the condition 
Wij = Wji has to be satisfied at least on average. In this 



case we can write 



Wn 



i) = 0, 



(12) 



and any initial conditions for the population size satisfies 
the stationary state. In the degree block variable rep- 
resentation we can recover the above condition by con- 
sidering a diffusion rate for each particle of the form 
Pk = Tk/Nk. The diffusion rate on any given edge from 
a subpopulation of degree /c to a subpopulation of degree 
k' is therefore given by 



dkk' = 



WQ{kk'f 



(13) 



and the degree block diffusion equations read in the un- 
correlated networks case as 



dtNkit) = -n + k'^^+''^ 



Wo 



(k) ■ 



(14) 



Since we know that by normalization Tf. = 
k^^^^^WQ{k^^^)/{k)^ we recover by definition the solution 
dtNk{t) = that allows any stationary value distribution 
N}^. Differently from the results obtained in the previous 
subsection, where each individual has the same probabil- 



ity p of leaving a subpopulation, eq. (12) shows that a 
population dependent diffusion process does not fix the 
subpopulation size, which can be given as a parameter 
of the model, with the only constraint that Nj^ > T/^, in 
order to make the diffusion process feasible. 



4. Epidemic spreading and the invasion threshold 

In order to explore the epidemic behavior in metapop- 
ulation models, the disease dynamics needs to be explic- 
itly considered inside each subpopulation. In the follow- 
ing we will consider the standard compartmentalization 
approach in which individuals exist in a certain num- 
ber of discrete states such as susceptible, infected or 
permanently recovered (Anderson & May 1992"). The 



paradigmatic epidemiological model one can consider, is 



the susceptible- infected-removed (SIR) model (Anderson 
fc May[ |1992[ |Murray[ |2QQ5[ ), where the total number of 
individuals Nj in the subpopulation j is partitioned in 
the compartment Sj{t)^ and Rj{t) denoting the 

number of susceptible, infected and recovered individ- 
uals at time t, respectively. By definition it follows 



Nj = Sj{t) + Ij{t) + Rj{t)' The disease transmission 
is described in an effective way. The probability that 
a susceptible individual acquires the infection from any 
given neighbor in an infinitesimal time interval dt is /3dt, 
where /3 defines the disease transmissibility. At the same 
time, infected vertices are cured and become recovered 
with probability p^dt. Individuals thus run stochastically 
through the susceptible — > infected recovered transi- 
tions, hence the name of the model. The SIR model as- 
sumes that recovered individuals are basically removed 
from the system, they do not participate anymore to the 
disease dynamics, due to their death or acquired immu- 
nization. Another popular model takes into account the 
possibility that infected individuals are again susceptible 
with probability jadt. In this case individuals thus run 
stochastically through the cycle susceptible infected 

susceptible, defining the so-called SIS model. The SIS 
model is mainly used as a paradigmatic model for the 
study of infectious diseases leading to an endemic state 
with a stationary and constant value for the prevalence of 
infected individuals, i.e. the degree to which the infection 
is widespread in the population. 

A basic parameter in the analysis of a single popula- 
tion epidemic outbreaks is the basic reproductive number 
Rq^ which counts the number of secondary infected cases 



generated by a primary infected individual (Anderson & 
May, 1992). Under the assumption of the homogeneous 
mixing of the population the basic reproductive number 
is defined as 

Ro = -. (15) 



It is straightforward to see from eq. (15) that in the sin- 
gle population case any epidemic will spread across a non 
zero fraction of the population only for Rq > 1. In this 
case the epidemic is able to generate a number of infected 
individuals larger than those who recover, leading to an 
increase in the overall number of infectious individuals 
I{t). The previous considerations lead to the definition of 
a crucial epidemiological concept - the epidemic thresh- 



old (Anderson & May, 1992). Indeed, if the spreading 



rate is not large enough to allow a reproductive number 
larger than one (i.e. /3 > /i), the epidemic outbreak will 
not affect a finite portion of the population and will die 
out in a finite amount of time. 

At the metapopulation level, however, the epidemic 
behavior on the global scale is determined also by the 
diffusion process of individuals. In particular, the effects 
due to the finite size of subpopulations and the stochastic 
nature of the diffusion might have a crucial role in the 
problem of resurgent epidemics, extinction and eradica- 



tion (Bah, 1997 Cross et a/., 


2005 Watts et a/., 


2005 


Vazquez 2007 


Cross et al 2007 


). Therefore it is impor- 



tant to consider the discrete nature of individuals. Indeed, 
each subpopulation may or may not transmit the infection 
to another subpopulation it is in contact with, depending 
on the occurrence or not of the travel event of at least one 
infected individual to the non-infected subpopulation dur- 
ing the entire epidemic evolution. The spreading process 
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across subpopulations will therefore occur with a prob- 
ability that is related to the diffusion probability of in- 
dividuals and the total number of individuals that will 



experience the infection (Ball 1997 Cross et 



2005 



2007|). In the case of epidemic processes with Rq < 1 
the epidemic will die with probability 1 and is not go- 
ing to spread across subpopulations. In a model like the 
SIS model, if Rq > I the number of infected individuals 
reaches a stationary state and the epidemic will eventually 
spread to different subpopulations, since locally endemic. 
In a model such as the SIR model, however, the epidemic 
within each subpopulation generates a finite fraction of 
infectious individuals in a given amount of time and even 
if Rq > 1 the diffusion rate must be large enough to en- 
sure the timely diffusion of infected individuals to other 
subpopulations of the metapopulation system, before the 
local epidemic outbreak dies out. This is captured by 
the definition of a new predictor of disease invasion, i?*, 
regulating the number of subpopulations that become in- 
fected from a single initially infected subpopulation; i.e. 
the analogous of the reproductive number at the subpop- 



ulation level (Bah, 1997 Cross et a/., 2005, 2007). 



This effect would not be captured by a simple de- 
terministic description that would allow any fraction pi 
of diffusing infected individual to inoculate the virus in 
a subpopulation not yet infected. In certain conditions 
this fraction pi may be a number smaller than one that 
just represents a mean- field average value. This is a 
common error of deterministic continuous approximations 
that allow the infection to persist and diffuse via "nano- 
individuals" that are not capturing the discrete nature 
of the real systems. For this reason, in the next section 
we will use an approach working at the level of subpop- 
ulations that allows to take into account effectively the 
fluctuations inherent to the diffusion process and the out- 
break extinction probability. 

4.1. GLOBAL INVASION THRESHOLD IN HOMOGENEOUS 
METAPOPULATION NETWORKS 

Let us consider a metapopulation system in which the 
initial condition is provided by a single introduction in 
a subpopulation of degree k and size A^^, given Rq > 1. 
While the stochastic nature of the process may lead in 
some cases to the extinction of the process, as Rq is above 
the epidemic threshold the epidemic will affect a finite 
fraction of the population with non zero probability. In 
the case of a macroscopic outbreak in a closed popula- 
tion the total number of infected individuals during the 
evolution of the epidemic will be equal to aN^ where a 
depends on the specific disease model used and the dis- 
ease parameter values. Each infected individual stays in 
the infectious state for an average time equal to the 
inverse of the recovery rate, during which it can travel 
to the neighboring subpopulation of degree with rate 
dkk'- To a first approximation we can therefore consider 
that the number of new seeds that may appear into a con- 
nected subpopulation of degree k' during the duration of 



the subpopulation epidemic is given by 



A 



kk' 



dkk' 



(16) 



In this perspective we can consider the metapopulation 
model in a coarse grained view (see Fig. |3| and provide a 
characterization of the invasion dynamics at the level of 
the subpopulations, translating epidemiological and de- 
mographic parameters into Levins-type metapopulation 
parameters of extinction and invasion rate. Let us define 
as the number of diseased subpopulation of degree 
k at generation 0, i.e. those which are experiencing an 
outbreak at the beginning of process. Each infected sub- 
population during the course of the outbreak will seed 
the infection in neighboring subpopulations defining the 
set D], of infected subpopulations at the following gener- 
ation and so on. This corresponds to a basic branching 



process (Harris, 1989 Bah, 1997 Vazquez, 2006) where 



the n— th generation of infected subpopulations of degree 
k is denoted DV:. 



diseased (i.e. infected) 
subopopulation 



metapopulation 




not infected subopopulation 



invasion dynamics at the 
subpopulation level: 
not infected subpop. O 
diseased subpop. ^ 



pop/ 



individuals level: 



O s 

of the % I 
• R 



Figure 3: Schematic representation of the invasion dynam- 
ics at the level of the subpopulations. The metapopulation 
system can be considered in a coarse grained perspective as 
a network where each node represents a subpopulation which 
can be infected (i.e. diseased) if it is reached by the virus as 
carried by the infected individuals diffusing on the system. 

In the early stage of the epidemics we assume that the 
number of subpopulations affected by an outbreak (with 
i?o > 1) is small and we can therefore study the evolu- 
tion of the number of diseased subpopulations by using 
a tree-like approximation relating D'^ with D^~^ . Let 
us first analyze the case of a metapopulation system in 
the form of a homogeneous random graph in which each 
subpopulation has the same degree k = k and population 
N. In this case we can drop the subscript index k (all 
subpopulations being equal) and obtain that 



D^'-^k-l) 



1 



1 

Ro 



1 



V 



(17) 



This equation assumes that each infected subpopulation 
of the (n — 1)— th generation, D^~^^ will seed with in- 
fected individuals a number of subpopulations during the 
course of the outbreak that depends on the product of: 
the number of neighbor subpopulations minus the one 
which originally transmitted the disease, ^ — 1, times the 
probability that the subpopulation is not already seeded 
by infected individuals, {1—D^~^)/V, and the probability 
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that the new seeded subp opulation wih 



break, i.e. ^1 — Rq ^^^^ (Bailey 
sion stems from the probabihty of 



experience an out- 
1975|. The last expres- 

: extinction P^xt = 1/^0 
given by a seed of a single infectious individual ( |Bailey[ 
|1975|) . The simplest case of homogeneous diffusion of in- 
dividuals dj^ = p/k yields A^^ = pNaji'^ /k. In order to 
obtain an explicit result we will consider in the following 
that i?o — 1 <^ 1, so that the system is poised very close to 
the epidemic threshold. In this limit we can approximate 
the outbreak probability as 



1 



^kki^o - 1), 



(18) 



and assuming that at the early stage of the epidemic 
D'^-^/V <C 1 we obtain 

D"" = pNa^-^^^^{Ro - l)D''-\ (19) 



This relation states that the number of subpopulations 
affected by an outbreak will increase only if the quantity 



R^ = pNa/jj 



-(i?o-l) > 1- 



(20) 



This relation defines the global invasion threshold^ i.e. the 
subpopulation reproductive number R^ that is the anal- 
ogous of the basic reproductive number Rq in structured 
metapopulation models. From the above expression it is 
possible to write the threshold condition on the mobility 
rate 



pN > 



k 



1 a 



1) 



(21) 



that fixes the threshold in the diffusion of individuals for 
the global spread of the epidemic in the metapopulation 
systems. In other words, this equation states that there 
is a minimum rate for the traveling of individuals in or- 
der to ensure that on average each subpopulation can 
seed more than one neighboring subpopulations. As it 



has been pointed out by Cross et al. (2007) we find that 



other factors such as the infectious period and the mobil- 
ity process are as much important as Rq in the spread of 
epidemics in structured populations. The constant a is 
larger than zero for any Rq > 1^ and in the SIR case for 
Rq close to 1 it can be approximated by ( [Murray , 2005): 



U., 2(J?o - 1) 

yielding the mobility threshold for the SIR model 

k iiRq 



pN > 



- 1 2(i?o - 1)2 ■ 



(22) 



(23) 



The above condition readily tells us that the closer to the 
epidemic threshold is the epidemic in the single subpop- 
ulation and the larger it has to be the traveling rate in 
order to sustain the global spread into the metapopulation 



model. We therefore find that we can define two differ- 
ent thresholds in a homogeneous metapopulation model. 
The first one is the local epidemic threshold Rq > 1 within 
each subpopulation and the second one R^ > 1 represents 
the global invasion threshold defining the traveling rate of 



individuals according to eq. (23). It is important to stress 
that when Rq increases the small Rq — ^ expansions are no 
longer valid and the invasion threshold is obtained only 
in the form of a complicate implicit expression. 

4.2. GLOBAL INVASION THRESHOLD IN 
METAPOPULATION NETWORKS WITH TRAFFIC 
DEPENDENT MOBILITY RATES 

The calculation for the global threshold becomes more 
complicate in the case of heterogeneous metapopulation 
networks. In this case eq. (17) has to include also the 



degree and population heterogeneities yielding: 

= ^ D]:r\k' - i)Xk'k{Ro - i)P{k\k') (i 



D 



n-l 



(24) 

This expression considers that each subpopulation of de- 
gree k' will seed the infection in a number k' — l of subpop- 
ulations corresponding to the number of neighboring sub- 
populations minus the one which originally transmitted 
the infection, the probability P{k\k') that each of the k' 
neighboring populations has degree /c, and the probability 
to observe an outbreak in the seeded population, where 
as before we considered the limit RQ — l<^lio obtain the 
the outbreak probability as Xk'k{Ro — !)• As in the pre- 
vious case of the homogeneous network, we consider the 
early stage of the epidemic in which (1 — D]^~^ /Vk) — 1. 
In addition we assume that degree correlations can be 
neglected and P{k\k') = kP{k)/{k) obtaining 



kP{k) 



{R^-i)Y,Dir\k'-i)\ 



k'k- 



(25) 



The behavior of the above expression depends on the spe- 
cific form of Xk'k that is determined by the diffusion rate 
dkk' ' Let us first consider the heterogeneous diffusion rate 
of Eq. ^ that gives 



A 



k'k 



p{k) a {kf 



k' 



p{k) 



{kk'fN, (26) 



where in the last expression we have considered that the 
size of each population of degree k is given by the station- 
ary diffusion process according to Eq. (10). The equa- 



tion describing the generation of infected subpopulations 
is therefore reading as 

Dk={Ro-l)^-^^^^^ (27) 



By defining 6"^ = D]^,k'^ {k' - 1) the last expression 
can be conveniently written in the iterative form 
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that allows the increasing of infected subpopulations and 
a global epidemic in the metapopulation process only if 

The subpopulation reproductive number is therefore an 
increasing function of the network heterogeneity that 
plays a role in the spread of the pathogen across sub- 
populations. In the case of an SIR epidemic within each 
subpopulation, the threshold on the mobility rate is pro- 
vided by the expression 



pN > 



(30) 



that differs from the homogeneous case for a correction 
factor depending on the topology of the network. Notice- 
ably, the ratio {k^+^)^ / - (k^^^^)) is extremely 
small in heavy-tailed networks and it is vanishing in the 
limit of infinite network size. This implies that the het- 
erogeneity of the metapopulation network is favoring the 
global spread of epidemics by lowering the global invasion 
threshold. 



4.3. GLOBAL INVASION THRESHOLD IN 
METAPOPULATION NETWORKS WITH POPULATION 
DEPENDENT MOBILITY RATES 

As a final case let us consider the realistic framework 
in which the diffusion rate of individuals is proportional 
to the ratio between traveling people and population size; 
i.e. pk = Tk/Nk- In subsection 3.2 we have seen that 
this case corresponds to the mean-field assumption in 
metapopulation models coupled by traveling fluxes, lead- 
ing to a stationary state in which the population Nj. is sta- 
tionary and independent on the diffusion process. Here 
\kk' = wo{kk')^ajj.~^ and by using the approximations 
considered in the previous cases the basic equation for 
the number of infected subpopulations reads as 



= (Ro - 1) 



k^+'^Pjk) 
{k) 



^D-,r'k^\k^-i). (31) 



Also in this case by using the auxiliary function O"^ = 
D'^,k'^(k' — 1) we obtain the recursive relation 

e- = (ii„-i)'^'"'' Hge--. (32) 

{k) /i 
yielding for the global invasion the condition 

^..W-i) <^"">:'^""' ^>L (as) 

{k) /i 

In the case of an SIR model for the intra-population dis- 
ease we obtain 



{k) 



where also in this case the mobility threshold is low- 
ered by the topological fluctuations of the network as the 
more heterogeneous is the metapopulation network and 
the smaller is the ratio {k)/{{P^'^^) - (fc^+^^)). In this 
respect it is worth remarking that in principle in an in- 
finite network with heavy-tails the mobility threshold is 
vanishing as the ratio (k) / {{k'^^'^^) - (A:^+^^)) 0. In an 
infinite network, however, the above equations should be 
rewritten in terms of the density of diseased subpopula- 
tions in order to avoid the pathological divergence of some 
terms. Finally, it is interesting to notice that the effect 
of the network heterogeneity on the subpopulation repro- 
ductive number is similar to that of the contact pattern 



heterogeneity on the basic reproductive number (Ander- 
son May] |1992| [Pastor-Satorras Vespignani| 1200101 
Lloyd fc May I |2001| [Barthelemy et aL\ |2005^ tressing 



even more the close analogy between the two metrics. 

4.4. LOCAL AND GLOBAL THRESHOLD IN REAL-WORLD 
CASES 

It is worth stressing that the previous expressions are 
approximate and valid only in the limit in which a small 
fraction of the populations in the system is affected and in 
which Rq — 1 <^ 1. It is however extremely relevant that 
metapopulations systems have intrinsically two epidemic 
thresholds. The emergence of a global epidemic is first 
constrained by the intrinsic epidemic threshold within 
each subpopulation, Rq > 1. If the epidemic process satis- 
fies this condition, each time an infectious individual seeds 
an epidemic within a subpopulation there is a finite prob- 
ability that a macroscopic fraction of the population will 
be affected by the outbreak. While this condition guaran- 
tees the intra-population spreading of the epidemic, the 
inter-population spreading is controlled by the coupling 
among subpopulations as quantified by the rate of diffus- 
ing/traveling individuals. The global invasion threshold 
condition R^ > 1 provides an estimate of the rate of dif- 
fusion of individuals above which the epidemic is able to 
affect a macroscopic fraction of the subpopulations defin- 
ing the meta-population network. 

At this point, it is useful to draw some gross estimate 
of the critical population coupling wq as a function of 
Ro and of a realistic value of /i. If we assume a very 
mild reproductive rate of about Rq :^ 1.1 and a value 
/i = 1/3 per unit time (1 day) as in many estimates for 
influenza strains, we obtain that wq per unit time must 
be larger than approximately 20 individuals per day in 
the case of homogeneous networks and one order of mag- 
nitude or more smaller in heterogeneous networks. This 
is a value met in most of the modern real transportation 
systems. For example, the world-wide air transportation 



network analyzed in ref. ( Barrat et al. , 2004| and briefly 
described in subsection 2.1 is characterized by a topology 
whose degree distribution moments which appear in the 
expressions of the global threshold are given by: {k) c^i 10 



(i^2+2^) - 2(i?o - 1)' 



(34) 



and - 7 • 10"^, given that 0.5 ( [Bar- 

rat et aL\ |2004[ ). Therefore the condition expressed in 
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eq. (34) states that an epidemic carried by air travelers 



would reach global proportion if the average number of 
travelers per day is larger than approximately 3 • 10~^, a 
constraint which is met in the airport network where the 
average daily traffic on a given connection has a minimum 
corresponding to :^ 10~^ The result of this estimate is 
in agreement with recent studies on contingency planning 
for a possible influenza pandemic, which show that travel 
restrictions, reducing the probability of any individual to 
leave an infected region, would not be able to consider- 
ably slow down the global spread unless > 90% or more 



effective ( Hollingsworth et a/., 2006 Cooper et a/., 2006 



Colizza et al] 2007 a[). Our analytical results shows that 



the invasion threshold is extremely small in realistic sit- 
uations and traffic reduction of more than one order of 
magnitude are in order to bring the system below the 
threshold. 



5. Epidemic behavior above the invasion 
threshold 

Above the invasion threshold we can assume that with 
finite probability the epidemics will affect a macroscopic 
fraction of subpopulations. In this limit, the stochastic 
effect due to the diffusion can be neglected and it is pos- 
sible to study the epidemic spreading in the system by 
the deterministic equations obtained from a mechanistic 
approach to the metapopulation model where the disease 
dynamics in each subpopulation can be viewed as a reac- 



tion process (Colizza et al. 20076). In the case of the SIR 



scheme the dynamics is identified by the following set of 
reaction equations: 



S 
I 



21 
R. 



(35) 
(36) 



In the SIS case the second reaction is just replaced by the 
reaction I ^ S. From the rate equations it is clear that 
the dynamics conserves the total number of individuals. 
Before the diffusion process, the Ij and Sj individuals be- 
longing to the same subpopulation j react according to 



the eqs. (35) and (36). In each node j the spontaneous 
process I ^ R simply consists in turning each Ij indi- 
vidual into an Rj individual with rate /i. This process 
account for the recovery of infected individuals from the 
disease. The process / + 5 ^ 2/ is related to the depen- 
dence of the transmissibility on the population density. In 
general, in large populations it is customary to consider 
that each individual has a finite number of contacts per 
unit time. In this case the probability that a suscepti- 
ble has a contact with an infectious individual is equal to 
the density of infectious individuals within the subpopu- 
lation j, i.e. Ij/Nj. If we consider a homogeneous mixing 
assumption within the population, the creation rate of in- 
fectious individuals will be provided by jSVj where Vj is 



an interaction kernel of the form 

r 



(37) 



It is natural also to consider different dependencies of the 
transmissibility with respect to the density, giving rise to 



different reaction kernels (Anderson & May, 1992 Col- 



izza et al., 20076). We will provide an analysis of the 



case of reaction kernels simulating population with inter- 
nal network structure and fully connected populations in 
a forthcoming paper (Colizza et al., in prep. 



5.1. DETERMINISTIC REACTION-DIFFUSION RATE 
EQUATIONS 

In order to provide the explicit equations describing 
the dynamical evolution of the metapopulation system, 
we generalize the basic mechanistic approach with degree 
block variables used in the previous section to the com- 
plete reaction diffusion process. We take into account 
the topological fluctuations of the coupling networks by 
introducing the quantities: 



Vk 



(38) 



which represent the the average number of / and S indi- 
viduals in subpopulations with degree k. Analogously, the 
reaction kernel in the homogeneous assumption is written 
for subpopulation in each degree block as = I^Sk/Nk. 
Again it is worth remarking that the degree block vari- 
ables assumes the statistical equivalence of subpopula- 
tions with the same degree k. While this approximation 
is in fair agreement with empirical analysis, real- world 
subpopulations have differences that the present analysis 
does not take into account. 

At the end of the reaction-diffusion process the varia- 
tion in the number of infected individuals in subpopula- 
tions of degree block k can be written as a discrete master 
equations with the form 



h{t^M)-h(t) 



Wr 



k ' 



(39) 



where the terms , identify the number of infected 
individuals that entered or left the class 1^ because of both 
the disease dynamics and the diffusion process. Assuming 
the general framework in which the diffusion probability 
out of a given subpopulation depends on its degree, p^, 
the depletion term can be simply evaluated as: 



wr 



Pkh + (1 -Pk)l~ih' 



(40) 



The depletion term is just the sum of the Ik individuals 
that diffuse away of the subpopulation (first term of the 
r.h.s.) and the infected individuals that do not diffuse 
away from the subpopulation but have a transition to the 
class Rk (second term of the r.h.s.). The positive term 
takes into account both the new infected individuals 
generated by the disease dynamics within the subpopu- 
lation and the infected individuals that diffuse from the 



■1-IATA, International Air Transport Association, http://www.iata.org/ 
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neighboring subpopulations with diffusion rate dkk'^ and 
it is given by: 



il-Pk)f3Tk + k }2 Pik'\k)dk'k [(1 - Ai)4' + PTk> 



(41) 

The first term on the r.h.s. considers the newly gener- 
ated infected individuals within the subpopulation with 
degree k. The factor I — Pk takes into account only those 
individuals that do not diffuse out of the subpopulation. 
The second term accounts for all the infected individuals 
arriving because of the diffusion process from neighbor- 
ing subpopulations. The factor k considers the presence 
of k neighboring subpopulations, each one contributing 
a fraction d^'k of its number of diffusing infected indi- 
viduals which is given by the non recovering plus the 
newly generated ones (1 — /i)//^/ + /3Tk'. Finally, on each 
connection edge we have to average over the probability 
that the neighboring subpopulation has degree A:', that is 
given by the weighted sum over the conditional probabil- 
ity P{k'\k). As for the simple diffusion in section 3, we 
consider an infinitesimal time interval At ^ and divide 



both terms of equation ( 39 ) to obtain the following set of 
differential equations 

dJk = -Pkh + (1 - Pk) [-/^h + P^k] + 

^kY,P{k'\k)dk^k [(1 - + pTk'] , (42) 

k' 

where all the parameter combinations are now infinitesi- 
mal transition rates. By considering the uncorrelated case 
P{k'\k) = k' P{k')/ (/c), we obtain the following dynamical 
reaction-rate equations 

dth = -Pkh + (1 - Pk) [-/J^h + P^k] + 

+ E ^'PWk'k [(1 - ^)h' + PTk'] . (43) 

Similar expressions can be written for the evolution of Sk 
and Rk as: 

dtSk = -PkSk - (1 - Pk)P^k + 

+ E ^'PWk'k [Sk' - pTk'] , (44) 



and 



dtRk = -PkRk + (1 - Pk)l~ih 



{k) 



Y,k'P{k')dk'kWk'^Rk']^ (45) 



5.2. THE EARLY STAGE OF THE EPIDEMIC OUTBREAK 

An explicit solution to the previous equations can be 
obtained for the early stages of the epidemic. In this case 
we can imagine to have a very small densities of infec- 
tious individuals in the metapopulation system so that in 



general we can neglect contributions of order In this 
setting the reaction kernel Vk can be approximated as 



{Nk — Ik — Rk)h 
Nk 



(46) 



where we have neglected all terms of order /| and consid- 
ered that Rk is of the same order of Ik in the early stage 
of the dynamics. The simplification of the reaction kernel 



allows to write eqs. (43) for the evolution of the density 



of infectious individuals in the following form: 



dth 



-Pkh 
k 



{l-pk){p-fi)Ik- 
Y,k'p{k')dk'k[(i-^- 



■p)ik'. 



(47) 



Explicit solutions to the above set of equations for the 
early dynamics of infectious individuals in degree block 
k can be found by considering the specific diffusion pro- 
cesses already introduced in the previous section. 

5.3. TRAFFIC DEPENDENT MOBILITY RATES 

By considering a uniform p and the expression for dkk' 
of eq. ^ we obtain after some simple algebra the follow- 
ing dynamical reaction-rate equations 



dth = -plk^{l-p){P-f^)lk^Pjj^^^ [(1 - /i + /3)7] , 

(48) 

that depend only on the densities of infectious individ- 
uals, and where I = Xl/c' ^(^O^fc'- ^ solution for the 
global density of infectious individual in the metapopu- 
lation system is obtained by averaging both terms of the 
equation over P(/c), obtaining: 



a,^p(fc)4 = aJ=(/3-/i)/, 



(49) 



where we have considered that P(k)k^^^^^ = {k^^^^^). 
This equation has the simple solution 



/ = /(0)e 



(/3-/i)^ 



(50) 



where 7(0) is the initial number of infected individuals 
in the metapopulation system. It readily states that the 
overall density of infectious individuals in the system can 
grow only if > /i, thus recovering the epidemic thresh- 
old condition Rq = jS / ji > 1. The metapopulation system 
exhibits at the deterministic level an epidemic threshold 
condition equivalent to that of each single population that 
sets the time scale for the whole system. Intuitively this 
is stating that if the epidemic is not able to proliferate in 
each local subpopulation, then it cannot produce a major 
outbreak at the metapopulation level. 

It is possible to solve the early time behavior for all 
subpopulations of a given degree block Ik{t) by plugging 
in the explicit solution of I{t) in eqs. (48). This yields: 



Ik{t) 



k^ 



,(/3-/i)^ ^^^g[(l-p)(/3-/i)- 



(51) 
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where A and Ck are parameters fixed by the initial con- 
ditions. If we assume that the metapopulation system is 
seeded with a total number of Iq initially infected individ- 
uals distributed homogeneously among subpopulations, 
i.e. 4(0) = lo/V = /(O) Vfc, we obtain: 

A = m and C, = m(^l-^^^ (52) 

If the Jo infected are distributed only in the fco— block 
subpopulations, /^(O) = (5fc,fcoJ(0)/P(fco), the coefficients 
A and Ck assume the following values: 



A = /(O) and C'k = /(O) 



kO 



Piko) 



(53) 



The above results show how the choice of initial condi- 
tions, whether if homogeneously distributed or locally dis- 
tributed, affects the early stages of the epidemic outbreak 
inside subpopulations of different block k. The change in 
the initial stage of the disease evolution in subpopula- 
tions depending on the degree block is confirmed by the 
numerical results reported in section 6. 

5.4. POPULATION DEPENDENT DIFFUSION RATE 

If we consider a population dependent diffusion rate, 
Pk = Tk/Nj., the system's behavior will be given by plug- 



ging into the set of eqs. (43) the degree dependent dif- 



fusion probability pk and the expression of the rate of 

diffusion on a link k' /c, dk'k = ^^n^'' ' -^^ ^P" 
proximation of early stage dynamics and considering the 
normalization condition Tk = k^~^^wo{k^~^^)/{k), we ob- 
tain: 



dth = -Pkh + (1 -Pk){P - lAh + 

^1+6' 



(54) 



where Vt = P{k)pklk' Proceeding along the line fol- 
lowed in the previous section, the solution can be found 
for the early stage behavior of the average number of infec- 
tious individuals I = P{^)^k by averaging both terms 
of the equation over P{k): 



yielding 



I = /(0)e 



i0-ii)t 



(55) 



(56) 



and thus recovering the epidemic threshold condition 
Rq = P/ ii> 1 also in this case. The early stage behavior 
of the epidemic in the deterministic approximation for the 
whole system does not differ from what observed in the 
previous heterogeneous frame. It is worth to mention that 
the solution for the population dependent diffusion rate 
is obtained under very general conditions for the size of 
the subpopulations N}. which can assume any value sub- 
ject only to the constraint N}. > Tj. to ensure a proper 
definition of the probability of diffusion. 



It is clear from the previous analysis that the deter- 
ministic equations are not capable to account for the inva- 
sion threshold as they consider the diffusion and reaction 
processes in a mean-field perspective that provides de- 
terministic equations for the average values. Indeed, the 
deterministic diffusion process allows any fraction pl^ of 
infected individual to deterministically seed new popula- 
tions, washing out the stochastic effects responsible for 
the invasion threshold. While these equations do not al- 
low to capture the invasion threshold, they provide a good 
description of the system and its behavior across degree 
classes above the invasion threshold, as we will show in 
the next section by comparing the analytical results with 
stochastic simulations at the mechanistic level. 



6. Mechanistic numerical simulations 

Here we provide extensive numerical simulations to 
support the theoretical picture described above. We re- 
port results from Monte Carlo simulations in a variety of 
different cases and compare them with the analytical find- 
ings. We adopt mechanistic numerical simulations where 
each single individual is tracked in time, during both the 
infection dynamics and the diffusion processes. The sys- 
tem evolves following a stochastic microscopic dynamics 
and at each time step it is possible to monitor quantities 
which depend on the subpopulation - such as e.g. the 
number Ij{t) of infectious individuals in the subpopula- 
tion j at time t - and also averages over blocks of nodes - 
e.g. the average number Ikif) of infected in subpopu- 
lations with degree k - or over the whole system, I{t). 
In addition, it is possible to study the evolution of the 
epidemic by monitoring the invasion dynamics at the lo- 
cal population level, and therefore measure the number 
of diseased subpopulations at time t, D{t). Given the 
stochastic nature of the dynamics, the experiment can be 
repeated with different realizations of the noise, different 
underlying graphs, and different initial conditions. This 
approach is equivalent to the real evolution of epidemic 
processes in the generated networks and can be used to 
validate the theoretical results obtained in the analytical 
approach. 

The substrate network is given by an uncorrelated 
complex network generated with the uncorrelated con- 



figuration model (Catanzaro et al 2005), based on the 



Molloy-Reed algorithm (Molloyfe Reed, 1995) with an ad- 



ditional constraint on the possible maximum value of the 
degree in order to avoid inherent structural correlations. 
The algorithm is defined as follows. Each node i is as- 
signed a degree ki obtained from a given degree sequence 
P{k) subject to the restriction ki < V^l^ . Here we as- 
sume a power-law degree distribution, P(k) ~ k~^ ^ with 
7 = 2.1 and 7 = 3. Links are then drawn to randomly 
connect pairs of nodes, respecting their degree and avoid- 
ing self-loops and multiple edges. Sizes oiV = 10^ and 
V = 10^ nodes have been considered. Weights on the con- 
nections among subpopulations are defined following the 
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statistical law found in real transportation system (see 
eq. ([T])). Therefore the weight on the link between sub- 
population i and subpopulation j is given by: 



Wi 



: wo{kikj 



(57) 



where ki and kj are the degrees of the subpopulations 
i and j, respectively. Here we fix k;o = 1, whereas 
assumes different values, including ^ = for uniform 
weights. This expression for the weights is then used to 
define the diffusion rates in the cases analytically investi- 
gated, as in eq. ^ and eq. (13). 

The dynamics proceeds in parallel and considers dis- 
crete time steps representing the unitary time scale r of 
the process. The reaction and diffusion rates are there- 
fore converted into probabilities and at each time step, 
the system is updated according to the following rules, 
a) Infection dynamics: i) The contagion process assumes 
that in each subpopulation j individuals homogeneously 
mix and have a finite number of contacts, so that the 
probability for a susceptible to contract a virus from an 
infected is proportional to the transmission rate and nor- 
malized to the subpopulation size, P/Nj. At each time 
step in the simulation, each susceptible is turned into an 
infectious with probability 1 — (1 — -^r)^^ . ii) At the same 
time, each infectious individual is subject to the recovery 
process and becomes recovered with probability /ir. b) 
After all nodes have been updated for the reaction, we 
simulate the diffusion process. Results shown in the fol- 
lowing subsections refer to the traffic dependent diffusion 
rate. 

6.1. GLOBAL AND LOCAL THRESHOLD IN 
HETEROGENEOUS METAPOPULATION MODELS 

In the previous sections we have shown that along with 
the usual local epidemic threshold Rq > 1^ the stochastic- 
ity and discreteness of the metapopulation diffusion pro- 
cess induce an intrinsic invasion threshold i?* > 1 - at 
the global level - which rules the invasion dynamics in 
the coarse-grained view of the system. This threshold de- 
termines whether the coupling between subpopulation is 
high enough in order to allow the virus to spread from 
one subpopulation to another and invade a finite portion 
of the whole system. Here we numerically investigate 
this phenomenon by studying a metapopulation model 
with heterogeneous structure {P{k) ^ k~^) and varying 
the coupling force between subpopulations. Initially, let 
us consider that the diffusion probability of an individ- 
ual on the heterogeneous metapopulation structure is lo- 
cally independent of the degree of the subpopulation - i.e. 
Pk = P - and heterogeneous on the links departing from 
a given subpopulation, following eq. The probability 
of diffusion from a subpopulation i to a subpopulation j 
for each individual in any given compartment located in 
i is therefore given by: 

wo{kikjy 



(58) 



at each time step, after the update for the local infec- 
tion dynamics within the subpopulations (see previous 
section), each individual in any compartment in subpop- 
ulation i moves to a neighboring subpopulation j with 
probability dij. We analyze different values of Rq by as- 
suming /i = 0.2 and /i = 0.02, and varying the value of the 
transmission rate (3. The simulations start with a local- 
ized initial condition given by the seeding of a subpopu- 
lation having degree /cq with /q = 10 infected individuals. 
This allows to monitor of the epidemic evolution in the 
metapopulation model and measure the final size of the 
epidemic, expressed in terms of the number or density of 
cases obtained in the whole system and the number of 
subpopulations experiencing an outbreak. 




300 



where = Wij represents the traffic in i. The sim- 
ulations proceeds according to the following procedure: 



300 



Figure 4: Metapopulation system's behavior above and be- 
low the global threshold. Results refer to i?o = 3, iV = 10^, 
and = 0.5. The epidemic is in both cases above the local 
threshold, leading to an exponential increase of /(t). Differ- 
ences in the diffusion probability values (top: p = 0.5, bottom: 
p = 10~^) show the effect of the global threshold on the num- 
ber D{t) of diseased subpopulations. D{t) is normalized to the 
system size V = 10^ for sake of visualization. 

In Figure [4] we analyze the behavior of I{t) and D(t) 
for = 3, above the local threshold, and for two values of 
the diffusion rate p = 0.5 and p = 10~^ that poise the sys- 
tem below and above the invasion threshold, respectively. 
While in both cases the figure shows an increase of the 
value of /(t), the behavior of D{t) is very different above 
and below the invasion threshold. Indeed, above the in- 
vasion threshold the number of affected subpopulations 
is increasing exponentially, while below the threshold the 
number of subpopulations remains small and goes to zero 
in a finite time. The increase in I{t) instead is guaran- 
teed also below the invasion threshold by the outbreak in 
the initial seeded population. On a longer time however 
I{t) keeps increasing only if the system is above the inva- 
sion threshold and new subpopulations are progressively 
infected. 
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Figure 5: Global threshold in a heterogeneous metapopulation system with traffic dependent diffusion rates. The left panel 
shows a 3D surface representing the value of the final epidemic size in the metapopulation system as a function of the local 
threshold Ro and of the diffusion probability p. If Ro approaches the threshold, larger values of the diffusion probability p need 
to be considered in order to observe a global outbreak in the metapopulation system. On the right, two plots showing the cross 
sections of the 3D plot at fixed values of Ro (top) and at fixed values of the traveling rate p (bottom) . 



While Figure [4] provides a clear evidence of the two 
separate threshold mechanisms, a complete analysis of 
the system phase diagram is obtained by analyzing the 
behavior of the global metapopulation attack rate, de- 
fined as the total fraction of cases R{oo)/N at the end of 
the epidemic, as a function of both Rq and p. In Figure [5) 
we report the global attack rate surface in the p-Ro space, 
and the two dimensional plots of the p and Rq crosscuts. 
Figure |5] clearly shows the effect of different couplings as 
expressed by the value of p in reducing the final size of 
the epidemic at a given fixed value of Rq . The smaller the 
value of Ro, the higher the coupling needs to be in order 
for the virus to successfully invade a finite fraction of the 
subpopulations, in agreement with the analytic result of 



eq. (30). This provides a clear illustration of the varying 



global invasion threshold as a function of the reproductive 
rate Rq. On the contrary, p-crosscuts show that whatever 
the value of p, Rq < 1 does not allow the epidemic to 
spread. 
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Figure 6: Effect of metapopulation structure heterogeneity 
on the global epidemic threshold. The final fraction of dis- 
eased subpopulations D(oo)/V at the end of the global epi- 
demic is shown as a function of the traveling diffusion rate 
p. A heterogeneous network with heavy-tailed degree distri- 
bution, P{k) ~ k~^'^ is compared to a homogeneous network 
with poissonian P{k) having the same size V — 10^ and same 



average degree. Here ^ = 0. 

Finally, it is possible to study the effect of the het- 
erogeneity of the metapopulation structure. Figure [6] 
shows the results obtained comparing a heterogeneous 
network characterized by a scale-free degree distribution 
P{k) ~ with a homogeneous network having the 

same size V = 10^ and same average degree. The pres- 
ence of topological fluctuations lead to a smaller ratio 
(A:i+^)V((A:2+2^) - {k^"^^^)), thus lowering the value of 
the mobility threshold with respect to the homogeneous 
network. 



6.2. EPIDEMICS ABOVE THE INVASION THRESHOLD 

Above the global invasion threshold R^ > 1, the epi- 
demic process is guaranteed to invade a macroscopic frac- 
tion of subpopulations and it is possible to inspect the 
validity of the results obtained in Section 4 with the de- 
terministic reaction-diffusion equations. A general con- 
clusion is that the global density of infectious individuals 
in the system in the early stage of the epidemic dynamics 
grows as 



m = me 



(/3-M)t 



(59) 



if the threshold condition, Rq = (3/ jj. > 1 is satisfied. The 
early time behavior expressed in the above equation is 
also independent of the parameters related to the diffusion 
process among subpopulations, such as the homogeneous 
diffusion rate p and the exponent which governs the re- 
lation between weights and subpopulation degrees. The 
analytic result of eq. ( 59 ) is confirmed in Fig. [7| where we 
show simulation results of the metapopulation epidemic 
model with traffic dependent diffusion rates. We con- 
sider systems of V = 10^ subpopulations each of initial 
size N = 10^, connected through a heterogeneous net- 
work having degree distribution P{k) k~^. The disease 
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parameters assume the following values: (3 = 0.04 and 
/i = 0.02, yielding Rq > 1. The simulations are seeded 
with 7(0) = 100 infectious individuals, homogeneously 
distributed among subpopulations. Both homogeneous 
{0 = 0) and heterogeneous {0 = 0.5) diffusions are consid- 
ered, as well as different values of the diffusion probability 
p = 0.5, 0.75, 1.0. 

Results in Fig. [7| show that the early behavior of the 
global density of infectious individuals is independent of 
the values of and p, and of the location of the initial 
seed, whether if homogeneously distributed among the 
subpopulations of a given degree block or localized in a 
single subpopulation. All simulations show an exponen- 
tial increase which confirms the analytic findings. 
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Figure 7: Evolution in time of the global density of infectious 
individuals in a heterogeneous metapopulation system with 
traffic dependent diffusion rates. Top: the effects of the hetero- 
geneity of diffusion (0), of the probability of diffusion (p), and 
of the distribution of initially infected individuals in the sys- 
tem (homogeneously distributed in a /c— degree block or con- 
centrated in a single subpopulation) are compared and found 
to produce the same early stage behavior. Bottom: changes 
in the initial condition value / produce the same exponential 
increase in the metapopulation system behavior. 

Numerical simulations also allow for the study of the 
dynamic behavior of infectious individuals in subpopula- 
tions of degree block k. The solutions obtained in subsec- 
tion 4.4 show a dependence of the early time behavior on 
the degree k of the subpopulation, pointing to a dynamics 
which switches on degree modes at different times. Re- 
sults of the numerical simulations confirm this findings, 
as shown in Fig. [Sj Here the disease parameters assume 
the same values as before, and the diffusion is governed by 
the values 6 = 0.5 and p = 0.75 for the numerical results 
reported in the top panel, whereas the effect of different 
values of is reported in the bottom panel. In order to see 
the effects of different initial conditions on the dynamical 



behavior of degree block k subpopulations (see eqs. (52) 
and (53)), we seed the epidemics with i) 10^ infectious 



individuals homogeneously distributed among subpopu- 
lations, or with ii) 10^ infectious individuals localized in 
subpopulations of degree block ko (results in the top panel 
correspond to ko = kmax)- While the global behavior 
I(t) is not affected by the choice of the initial conditions 
(see previous Figure), the subpopulations experience out- 
breaks at different times, as brought and delayed by the 
diffusion dynamics. The system heterogeneity, as con- 
tained in the factor k^^^ / {k^^^) of the explicit solution of 
//e(t), differentiates the evolution of the degree block sub- 
populations at short times. Numerical results obtained 
for the study of the effect of traffic heterogeneity (Fig. [s] 
bottom) are compared with the analytical findings of sub- 
section 4.4. 
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Figure 8: Evolution in time of the density of infectious indi- 
viduals in degree block k subpopulations. Top: changes in the 
initial conditions (homogeneous vs. localized in kmax) yield 
different behaviors in the early time dynamics of distinct de- 
gree blocks. Bottom: changes in the value of impact differ- 
ently the time evolution of the average number of infected Ik 
in each degree block, whether if /c = /co (i.e. kmin) or /c / /co 
(i.e. k — kmax)- 



7. Conclusions and outlook 

Here we have introduced an analytic framework in 
terms of degree block variables which allows to gain in- 
sights on the behavior of mechanistic metapopulation epi- 
demic models which explicitly include demographic and 
mobility heterogeneities. The system is shown to display 
a local epidemic threshold which depends on the disease 
parameter values only and is responsible for the epidemic 
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outbreak at the local scale, and a global epidemic thresh- 
old which determines the invasion dynamics of the sub- 
populations and depends critically on the disease param- 
eters and the diffusion rates of the individuals. Changes 
in coupling between subpopulations are shown to have 
critical implications for disease extinction. 

The results provide useful insights for the basic theo- 
retical understanding of mechanistic epidemic models in 
complex environments, which can then be used to build 
more realistic data-driven large-scale computational ap- 
proaches for real case scenarios and spatially targeted 
control measures. However, several key theoretical and 
practical issues are still to be addressed. Data on human 
dynamics at the local level, i.e. within any subpopula- 
tion, could push forward a more sophisticate theoretical 
framework for the local infection dynamics, to go beyond 



the homogeneous mixing assumption ( [Meyers et al , 2005 



Lloyd-Smith et al. 2005). The behavior of metapopula- 



tion models characterized by complex internal structure 
in each patch is a major question for the theoretical epi- 
demiology of the future. In addition, more realistic and 
detailed diffusion patterns should be included in order to 
better model the coupling terms by including non-Markov 
processes and introducing elements of memory in the sys- 
tem. Obviously this corresponds to the need for more 
accurate data on population behavior, such as fraction 
of commuters, probability of short /medium/long range 
travel, trip duration, and so on ( [Riley [ [2 00 7[ ). Additional 
levels of heterogeneity can be also included in the diffusive 
patterns by introducing a dependence of the probability 
of diffusion on the stage of the disease. In many real cases, 
e.g. the severity of symptoms or hospitalization measures 
would prevent the diffusion out of a patch to a portion 
of the population. The impact of heterogeneities in trav- 
eling pattern of individuals depending on their infection 
state could provide additional insights fundamental to the 
study of global extinction and eradication. All these im- 
provements and future directions would help filling the 
gap between the evidence from increasingly realistic epi- 
demic models and their theoretical understanding. 
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